*** to be run after 1_Table_2\7_estimation_sample_and_variables.do

cap log close	

	
log  using  .../AA_8_trim_07Sep2020.smcl,replace

version 13	
use .../sample_estimation_1_03Sept2020.dta, clear



preserve
tempfile tmp
expand 10 if pweight == 10
sum netincdiff if sedum == 1, detail
gen k = 0 
replace k = (netincdiff<r(p1)  |netincdiff>r(p99)) if sedum == 1
bysort id: egen kk = max(k)

sum netincdiff if sedum == 0, detail
gen k1 = 0
replace k1 = (netincdiff<r(p1)  |netincdiff>r(p99)) if sedum == 0
bysort id: egen kk1 = max(k1)

by id: keep if _n==1
keep if kk==1 | kk1 == 1
keep id kk1 kk 
sort id 
save `tmp'
restore 
sort id year
merge m:1 id  using `tmp'
drop if kk == 1 | kk1== 1
drop kk kk1


*** 1 ***	
	cap program drop prog2_lnf
	set type double					

	program define prog2_lnf

	args todo b lnf

	tempvar ISE IU ISELF IULF logitp2 logitp3 logitp4 b1 b2 b3 b4   a1 a2 a3 a4 logitp5 logitp6 logitp7 logitp8  lnfie last lE1 lE2 lE3   lE1LF lE2LF lE2LF ///
	        lU1  lU2 lU3 lU1LF lU2LF lU3LF l1 ll1 l2 ll2 l3 ll3 l4 ll4 l5 ll5 l6 ll6 l7 ll7 l8 ll8  mlSE1 mlSE2  mlSE3  mlE1 mlE2 mlE3 mlSELF1  mlSELF2 mlSELF3  mlELF1 mlELF2 mlELF3 ///
	        h1SE h2SE h3SE  sum1SE sum2SE  sum3SE  ST1SE ST2SE ST3SE lE3LF h1SELF h2SELF  h3SELF  sum1SELF sum2SELF  sum3SELF  ST1SELF ST2SELF ST3SELF ///
		h1U h2U  h3U  sum1U sum2U  sum3U  ST1U ST2U ST3U h1ULF h2ULF   h3ULF  sum1ULF sum2ULF  sum3ULF  ST1ULF ST2ULF ST3ULF lastSE lastU lastSELF lastULF lnfi
		
	tempname p2 p3 p4 p1
	
        *defining equations
    mleval `ISE'     = `b', eq(1)
	mleval `IU'      = `b', eq(2)
    mleval `ISELF'   = `b', eq(3)
	mleval `IULF'    = `b', eq(4)
	mleval `a1'       = `b', scalar eq(5)
	mleval `a2'       = `b', scalar eq(6)
	mleval `b1'       = `b', scalar eq(7)
	mleval `b2'       = `b', scalar eq(8)
	mleval `logitp2' = `b', scalar  eq(9) //5
*	mleval `logitp3' = `b', scalar  eq(10) //5



	

	quietly {
			* to make sure probabilities constrained between 0 and 1
			* for identification, p11==1-p21-p22 and m11==0
			scalar `p2' = exp(`logitp2')/(1+exp(`logitp2'))
			scalar `p1' = 1             /(1+exp(`logitp2'))

			***** for SELFEMPLOYED spells
			
			* hazard functions for the three types
			by $S_E_id $spell: gen double `h1SE' = 1-exp(-exp(`ISE' + `a1'))	if $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `h2SE' = 1-exp(-exp(`ISE' ))			if $SE_E ==1 & $LF == 0
			
                        * first part of survival function 
			by $S_E_id $spell: gen double `sum1SE' = sum(exp(`ISE' + `a1'))		if $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `sum2SE' = sum(exp(`ISE' ))			if $SE_E ==1 & $LF == 0
		
			
			* cumulative product of survival functions
			by $S_E_id $spell: gen double `ST1SE' = exp(-`sum1SE'[_N])			if _n==_N & $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `ST2SE' = exp(-`sum2SE'[_N])			if _n==_N & $SE_E ==1 & $LF == 0
			

			* identify last obs of the spell
			by $S_E_id $spell: gen byte `lastSE' = (_n==_N) if $SE_E ==1 & $LF == 0
			replace `lastSE'=0 if `lastSE'==.

			* likelihood functions of the spells for each type
			by $S_E_id $spell: gen double `lE1' =  `ST1SE' * ((`h1SE'/(1-`h1SE'))^$event ) if _n==_N & $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `lE2' =  `ST2SE' * ((`h2SE'/(1-`h2SE'))^$event ) if _n==_N & $SE_E ==1 & $LF == 0
			
			
			replace `lE1'=0 if `lE1'==.
			replace `lE2'=0 if `lE2'==.
			

			******** for EMPLOYED 
			by $S_E_id $spell: gen double `h1U' = 1-exp(-exp(`IU' + `b1' )) if $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `h2U' = 1-exp(-exp(`IU' )) if $SE_E ==0 & $LF == 0
			
			
			by $S_E_id $spell: gen double `sum1U' = sum(exp(`IU' + `b1' )) if $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `sum2U' = sum(exp(`IU' )) if $SE_E ==0 & $LF == 0
			

			by $S_E_id $spell: gen double `ST1U' = exp(-`sum1U'[_N]) if _n==_N & $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `ST2U' = exp(-`sum2U'[_N]) if _n==_N & $SE_E ==0 & $LF == 0
			
			
			by $S_E_id $spell: gen byte `lastU' = (_n==_N) if $SE_E ==0 & $LF == 0
			replace `lastU'=0 if `lastU'==.

			by $S_E_id $spell: gen double `lU1' =  `ST1U' * ((`h1U'/(1-`h1U'))^$event ) if _n==_N & $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `lU2' =  `ST2U' * ((`h2U'/(1-`h2U'))^$event ) if _n==_N & $SE_E ==0 & $LF == 0
                        
			               
			replace `lU1'=0 if `lU1'==.
			replace `lU2'=0 if `lU2'==.
			
			
			***** for SELFEMPLOYED left-censored spells
			
			* hazard functions for the three types
			by $S_E_id $spell: gen double `h1SELF' = 1-exp(-exp(`ISELF' + `a2'))   if $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `h2SELF' = 1-exp(-exp(`ISELF'))   if $SE_E ==1 & $LF == 1
			
			
                        * first part of survival function
			by $S_E_id $spell: gen double `sum1SELF' = sum(exp(`ISELF' + `a2'))  if $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `sum2SELF' = sum(exp(`ISELF' ))    if $SE_E ==1 & $LF == 1
			
			
			* cumulative product of survival functions
			by $S_E_id $spell: gen double `ST1SELF' = exp(-`sum1SELF'[_N]) if _n==_N & $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `ST2SELF' = exp(-`sum2SELF'[_N]) if _n==_N & $SE_E ==1 & $LF == 1
			
			
			* identify last obs of the spell
			by $S_E_id $spell: gen byte `lastSELF' = (_n==_N) if $SE_E ==1 & $LF == 1
			replace `lastSELF'=0 if `lastSELF'==.

			* likelihood functions of the spells for each type
			by $S_E_id $spell: gen double `lE1LF' =  `ST1SELF' * ((`h1SELF'/(1-`h1SELF'))^$event ) if _n==_N & $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `lE2LF' =  `ST2SELF' * ((`h2SELF'/(1-`h2SELF'))^$event ) if _n==_N & $SE_E ==1 & $LF == 1
			
           
			replace `lE1LF'=0 if `lE1LF'==.
			replace `lE2LF'=0 if `lE2LF'==.
			
			

			******** for EMPLOYED left-censored spells
			by $S_E_id $spell: gen double `h1ULF' = 1-exp(-exp(`IULF' + `b2')) if $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `h2ULF' = 1-exp(-exp(`IULF' )) if $SE_E ==0 & $LF == 1
			

			by $S_E_id $spell: gen double `sum1ULF' = sum(exp(`IULF' + `b2')) if $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `sum2ULF' = sum(exp(`IULF' )) if $SE_E ==0 & $LF == 1
			
			
			by $S_E_id $spell: gen double `ST1ULF' = exp(-`sum1ULF'[_N]) if _n==_N & $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `ST2ULF' = exp(-`sum2ULF'[_N]) if _n==_N & $SE_E ==0 & $LF == 1
			
			
			by $S_E_id $spell: gen byte `lastULF' = (_n==_N) if $SE_E ==0 & $LF == 1
			replace `lastULF'=0 if `lastULF'==.

			by $S_E_id $spell: gen double `lU1LF' =  `ST1ULF' * ((`h1ULF'/(1-`h1ULF'))^$event ) if _n==_N & $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `lU2LF' =  `ST2ULF' * ((`h2ULF'/(1-`h2ULF'))^$event ) if _n==_N & $SE_E ==0 & $LF == 1
			
                                       
			replace `lU1LF'=0 if `lU1LF'==.
			replace `lU2LF'=0 if `lU2LF'==.
			
	    
				// individual log likelyhood for each type
			by $S_E_id: gen double `l1'=`lE1' +`lU1'+`lU1LF'+`lE1LF'
			by $S_E_id: gen double `ll1'=`p1'*exp(sum(log(`l1')))
			by $S_E_id: gen double `l2'=`lU2'+`lE2' +`lU2LF'+`lE2LF'
			by $S_E_id: gen double `ll2'=`p2'*exp(sum(log(`l2')))




			by $S_E_id: gen byte `last' = (_n==_N) 
	
			* individual contribution to the log likelihood 
			gen double `lnfie' = cond( !`last'  , 0, $weight * ln(( `ll1') + ( `ll2')   ))
				
					
			mlsum `lnf' = `lnfie'

	}
end

	
	
	

set type double		


tempvar mysamp 
		tempname b b01 b02 b1 b2 b03 b04 b3 b4 lnf V  a1 a2 a3  a4 b4 masshE masslE masshSE masslSE ///
		       logitp2 logitp3 logitp4  masslSELF masslELF  masshSELF masshELF

  *now put all the variables togetherlocal
 *   se_netincdiff se_convexity se_female_1 se_ageb se_edu_2-se_edu_3  se_yr_2-se_yr_4 se_reg_2-se_reg_5
 *    lnt_se t_se t2_se 
 *    e_netincdiff e_convexity e_female_1 e_ageb e_edu_2-e_edu_3  e_yr_2-e_yr_4 e_reg_2-e_reg_5
 *  lnt_e  t_e t2_e
 *  self_netincdiff self_convexity  self_female_1 self_ageb self_edu_2-self_edu_3 self_yr_2-self_yr_4 self_reg_2-self_reg_5
 * t_self t2_self lnt_self
 * elf_netincdiff elf_convexity elf_female_1 elf_ageb elf_edu_2-elf_edu_3  elf_yr_2-elf_yr_4 elf_reg_2-elf_reg_5
 * t_elf t2_elf lnt_elf
global varse " se_netincdiff se_convexity se_female_1 se_ageb se_edu_2-se_edu_3  se_yr_2-se_yr_4 se_reg_2-se_reg_5"
global vare  " e_netincdiff e_convexity e_female_1 e_ageb e_edu_2-e_edu_3  e_yr_2-e_yr_4 e_reg_2-e_reg_5"
global varself " self_netincdiff self_convexity  self_female_1 self_ageb self_edu_2-self_edu_3 self_yr_2-self_yr_4 self_reg_2-self_reg_5"
global varelf  "elf_netincdiff elf_convexity elf_female_1 elf_ageb elf_edu_2-elf_edu_3  elf_yr_2-elf_yr_4 elf_reg_2-elf_reg_5"
global dumse "  lnt_se"
global dume "    lnt_e"
global dumself " lnt_self "
global dumelf " lnt_elf"

 
set more off	
// initial values for ll function taken from clogclog		
glm durvar  $varse  $dumse if sedum==1 & lcens_spell ==0 ,  f(b) l(cloglog) 
matrix `b01' = e(b)
matrix coleq `b01' = hazardE
	
glm durvar   $vare  $dume  if sedum==0 & lcens_spell ==0   ,  f(b) l(cloglog) 	
matrix `b02' = e(b)
matrix coleq `b02' = hazardU

glm durvar  $varself  $dumself  if sedum==1 & lcens_spell ==1 ,  f(b) l(cloglog) 
matrix `b03' = e(b)

matrix coleq `b03' = hazardELF
	
glm durvar   $varelf  $dumelf  if sedum==0 & lcens_spell ==1 ,  f(b) l(cloglog) 	
matrix `b04' = e(b)

matrix coleq `b04' = hazardULF

local LL1 = e(ll)
// create matrix with initial betas from cloglog and initial mass point location and probability (given)

matrix `a1' = (-0.1 )
matrix colnames `a1' = a1:_cons

matrix `a2' = (-1.3)
matrix colnames `a2' = a2:_cons

matrix `a3' = (-3.2)
matrix colnames `a3' = a3:_cons

matrix `b1' = (-0.2 )
matrix colnames `b1' = b1:_cons

matrix `b2' = (-1.2 )
matrix colnames `b2' = b2:_cons
matrix `b3' = (-1.5 )
matrix colnames `b3' = b3:_cons


matrix `logitp2' = (0.5)
matrix colnames `logitp2' = logitp2:_cons
*matrix `logitp3' = (-2.1)
*matrix colnames `logitp3' = logitp3:_cons
*matrix `logitp4' = (-1.1)
*matrix colnames `logitp4' = logitp4:_cons


matrix `b1' = `b01',`b02',`b03',`b04',`a1',`a2',`b1',`b2',`logitp2'


matrix list `b1'

// define macro that identify id variable
global S_E_id "id"
// define macro that identify se or E spell
global SE_E "sedum"
// define macro that identy the spell
global spell "a_spell"
// define macro that identify event
global event "durvar" 
 // define macro that identify left_censored spells
global LF "lcens_spell" 
 
// define weight
global weight "pweight"
 
 
sort id a_spell year
 
*cap log close
*log  using /ssb/ovibos/a2/ekstern/apa/prog/apa/survival/moredefition/durmod_`gruppo'.log,replace

//ml model, d0, no gradient or Hessian provided
ml model d0 prog2_lnf (hazardE: durvar = $varse  $dumse ) ///
                       (hazardU: durvar =  $vare  $dume  ) ///
		       (hazardELF: durvar = $varself $dumself   ) ///
		       (hazardULF: durvar = $varelf $dumelf  ) ///
					    (a1: ) (a2: ) (b1: ) (b2: )    ///
					    (logitp2: )   ///
                                , waldtest(0) //init(`b1')   

ml init hazardE:se_netincdiff    		=  -.3329825 
ml init hazardE:se_convexity    		=  .064971   
ml init hazardE:se_female_1    			=  -.044343  
ml init hazardE:se_ageb    				=  -.013019  
ml init hazardE:se_edu_2    			=  .018278   
ml init hazardE:se_edu_3    			=  .167996   
ml init hazardE:se_yr_2    				=  -.0840794 
ml init hazardE:se_yr_3    				=  .0176366 
ml init hazardE:se_yr_4    				=  -.2360357 
ml init hazardE:se_reg_2    			=  .0371971  
ml init hazardE:se_reg_3    			=  .021455   
ml init hazardE:se_reg_4    			=  .0504551  
ml init hazardE:se_reg_5    			=  .0962647  
ml init hazardE:lnt_se    				=  -.5238554 
ml init hazardE:_cons    				=  -1.023804 
ml init hazardU:e_netincdiff    		=  2.273436 
ml init hazardU:e_convexity    			=  -.2202451 
ml init hazardU:e_female_1    			=  .6399796  
ml init hazardU:e_ageb    				=  .0270455  
ml init hazardU:e_edu_2    				=  .1575     
ml init hazardU:e_edu_3    				=  .0366636  
ml init hazardU:e_yr_2    				=  .2521497  
ml init hazardU:e_yr_3    				=  .3393875  
ml init hazardU:e_yr_4    				=  -.1081969 
ml init hazardU:e_reg_2    				=  -.2560685 
ml init hazardU:e_reg_3    				=  -.3333598 
ml init hazardU:e_reg_4    				=  -.4136283 
ml init hazardU:e_reg_5    				=  -.4635098 
ml init hazardU:lnt_e    				=  -.4962889 
ml init hazardU:_cons    				=  -2.823308 
ml init hazardELF:self_netincdiff    	=  -.8615843 
ml init hazardELF:self_convexity    	=  -.0607664 
ml init hazardELF:self_female_1    		=  .2182683  
ml init hazardELF:self_ageb    			=  -.0358902 
ml init hazardELF:self_edu_2    		=  -.1504745 
ml init hazardELF:self_edu_3    		=  .0290571  
ml init hazardELF:self_yr_2    			=  -.1833268 
ml init hazardELF:self_yr_3    			=  -.3586449 
ml init hazardELF:self_yr_4    			=  -.79577   
ml init hazardELF:self_reg_2    		=  .0344597 
ml init hazardELF:self_reg_3    		=  -.0563954 
ml init hazardELF:self_reg_4    		=  .1703456  
ml init hazardELF:self_reg_5    		=  .1289736  
ml init hazardELF:lnt_self    			=  .010052   
ml init hazardELF:_cons    				=  -.884399  
ml init hazardULF:elf_netincdiff    	=  3.005003  
ml init hazardULF:elf_convexity    		=  -.2365221 
ml init hazardULF:elf_female_1    		=  .8126965 
ml init hazardULF:elf_ageb    			=  -.0456921 
ml init hazardULF:elf_edu_2    			=  .079268   
ml init hazardULF:elf_edu_3    			=  .217938   
ml init hazardULF:elf_yr_2    			=  .1006784  
ml init hazardULF:elf_yr_3    			=  -.0234075 
ml init hazardULF:elf_yr_4    			=  -.7619868 
ml init hazardULF:elf_reg_2    			=  -.4285836 
ml init hazardULF:elf_reg_3    			=  -.5198312 
ml init hazardULF:elf_reg_4    			=  -.5543597 
ml init hazardULF:elf_reg_5    			=  -.5843642 
ml init hazardULF:lnt_elf    			=  -.1714959 
ml init hazardULF:_cons    				=  -1.53394  
ml init a1:_cons    					=  -.5069521
ml init a2:_cons    					=  -1.328561  
ml init b1:_cons    					=  -2.774351 
ml init b2:_cons    					=  -1.839358 
ml init logitp2:_cons    				=  -1.522019
					
ml maximize, trace  search(off) difficult gradient

 
 eststo m

local gruppo = 1

eststo r1lr: nlcom (p1: 1/(1+exp([logitp2]_cons))) ///
		   (p2: exp([logitp2]_cons)/(1+exp([logitp2]_cons))), post

		   
log close


